Advanced statistics for physics analysis - Exercise 8

Exercise 1

Students from the Bachelor degree in Physics performed an experiment to study the Zeeman effect. The apparatus contains a Ne source lamp whose position can be changed. During the setting up of the apparatus, the source position has to be adjusted in order to maximize the intensity of the detected light signal. The following table gives the position of the source (in mm) and the corresponding height of the peak (arbitrary units) for the wavelength under study:

\(x_i\) = 2.44, 3.49, 3.78, 3.31, 3.18, 3.15, 3.1, 3.0, 3.6, 3.4

\(y_i\) = 129, 464, 189, 562, 589, 598, 606, 562, 360, 494

Assume a quadratic dependence of the peak height, \(y_i\) as a function of the source position \(x_i\),

\[ f(x)=c_0 + c_1x + c_2x^2 \]

All the measured values are affected by a Gaussian noise with zero mean, such that

\[ y_i = f(x_i) + \epsilon \]

where \(\epsilon\) follows a normal distribution with mean \(\mu = 0\) and unknown standard deviation, \(\sigma\).

1. Build a Markov Chain Monte Carlo to estimate the best parameters of the quadratic dependence of the data and the noise that affects the measured data.

# This code is a modified version of the code that can be found at https://github.com/ehalley/PBI/tree/master/PBI_scripts

library(mvtnorm)

# Metropolis (MCMC) algorithm to sample from function func.
# The first argument of func must be a real vector of parameters, 
# the initial values of which are provided by the real vector thetaInit.
# func() returns a two-element vector, the logPrior and logLike 
# (log base 10), the sum of which is taken to be the log of the density 
# function (i.e. unnormalized posterior). The MCMC sampling PDF is the 
# multivariate Gaussian with fixed covariance, sampleCov. A total of 
# Nburnin+Nsamp samples are drawn, of which the last Nsamp are kept. 
# ... is used to pass data, prior parameters etc. to func().
# If demo=FALSE (default), then
# return a Nsamp * (2+Ntheta) matrix (no names), where the columns are
# 1:  log10 prior PDF
# 2:  log10 likelihood
# 3+: Ntheta parameters
# (The order of the parameters in thetaInit and sampleCov must match.)
# If demo=TRUE, return the above (funcSamp) as well as thetaPropAll, a 
# Nsamp * Ntheta matrix of proposed steps, as a two element named list.
metrop <- function(func, thetaInit, Nburnin, Nsamp, sampleCov, verbose, 
                   demo=FALSE, ...) {

  Ntheta   <- length(thetaInit)
  thetaCur <- thetaInit
  funcCur  <- func(thetaInit, ...) # log10
  funcSamp <- matrix(data=NA, nrow=Nsamp, ncol=2+Ntheta) 
  # funcSamp will be filled and returned
  nAccept  <- 0
  acceptRate <- 0
  
  if(demo) {
    thetaPropAll <- matrix(data=NA, nrow=Nsamp, ncol=Ntheta)
  }
  
  for(n in 1:(Nburnin+Nsamp)) {

    #cat("cur", funcCur)
    # Metropolis algorithm. No Hastings factor for symmetric proposal
    if(is.null(dim(sampleCov))) { # theta and sampleCov are scalars
      thetaProp <- rnorm(n=1, mean=thetaCur, sd=sqrt(sampleCov))
    } else {
      thetaProp <- rmvnorm(n=1, mean=thetaCur, sigma=sampleCov, 
                           method="eigen")
    }
    funcProp  <- func(thetaProp, ...) 
    logMR <- sum(funcProp) - sum(funcCur) # log10 of the Metropolis ratio
    if(logMR>=0 || logMR>log10(runif(1, min=0, max=1))) {
      thetaCur   <- thetaProp
      funcCur    <- funcProp
      nAccept    <- nAccept + 1
      acceptRate <- nAccept/n
    }
    if(n>Nburnin) {
      funcSamp[n-Nburnin,1:2] <- funcCur
      funcSamp[n-Nburnin,3:(2+Ntheta)] <- thetaCur
      if(demo) {
        thetaPropAll[n-Nburnin,1:Ntheta] <- thetaProp
      }
    }
    if( is.finite(verbose) && (n%%verbose==0 || n==Nburnin+Nsamp) ) {
      s1 <- noquote(formatC(n,          format="d", digits=5, flag=""))
      s2 <- noquote(formatC(Nburnin,    format="g", digits=5, flag=""))
      s3 <- noquote(formatC(Nsamp,      format="g", digits=5, flag=""))
      s4 <- noquote(formatC(acceptRate, format="f", digits=4, width=7, 
                            flag=""))
      cat(s1, "of", s2, "+", s3, s4, "\n")
    }
  }
  if(demo) {
    return(list(funcSamp=funcSamp, thetaPropAll=thetaPropAll))
  } else {
    return(funcSamp)
  }
}

##### Functions to provide evaluations of prior, likelihood and posterior for the
##### quadratic model, plus sampling from the prior

# theta is vector of parameters; obsdata is 2 column dataframe with names [x,y].

# Return c(log10(prior), log10(likelihood)) (each generally unnormalized) of the quadratic model
logpost.quadraticmodel <- function(theta, obsdata) {
  logprior <- logprior.quadraticmodel(theta)
  if(is.finite(logprior)) { 
    return( c(logprior, loglike.quadraticmodel(theta, obsdata)) )
  } else {
    cat("logprior is not finite!")
    return( c(-Inf, -Inf) )
  }
}

# Return log10(likelihood) for parameters theta and obsdata
# dnorm(..., log=TRUE) returns log base e, so multiply by 1/ln(10) = 0.4342945
# to get log base 10
loglike.quadraticmodel <- function(theta, obsdata) {
  # convert alpha to b_1 and log10(ysig) to ysig
  theta[2] <- tan(theta[2])
  theta[4] <- 10^theta[4]
  modPred <- drop( theta[1:3] %*% t(cbind(1,obsdata$x,obsdata$x^2)) )
  # Dimensions in above mixed vector/matrix multiplication: [Ndat] = [P] %*% [P x Ndat] 
  logLike <- (1/log(10))*sum( dnorm(modPred - obsdata$y, mean=0, sd=theta[4], log=TRUE) )
  return(logLike)
}

# Return log10(unnormalized prior)
logprior.quadraticmodel <- function(theta) {
  b0Prior      <- dnorm(theta[1], mean=0, sd=10)
  alphaPrior   <- 1
  b2Prior      <- dnorm(theta[3], mean=0, sd=5)
  logysigPrior <- 1 
  logPrior <- sum( log10(b0Prior), log10(alphaPrior), log10(b2Prior), log10(logysigPrior) )
  return(logPrior)
}

library(gplots) 

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

### Define true model and simulate experimental data from it

set.seed(57)
Ndat <- 20
xrange <- c(0,10)
sigTrue <- 2
modMat <- c(0.69295, -0.48641, -0.76994)

x.in <- c(2.44, 3.49, 3.78, 3.31, 3.18, 3.15, 3.1, 3.0, 3.6, 3.4)
y.in <- c(129, 464, 189, 562, 589, 598, 606, 562, 360, 494)
x <- (x.in - mean(x.in))/sd(x.in)
y <- (y.in - mean(y.in))/sd(y.in)

plot(x, y,
     col = "red3",
     main = "Zeeman data")

# True parameters, transformed to be conformable with model to be used below
thetaTrue <- c(modMat[1], atan(modMat[2]), modMat[3], log10(sigTrue))
obsdata <- data.frame(cbind(x,y)) # columns must be named "x" and "y"

### Define model and infer the posterior PDF over its parameters

# Model to infer: linear regression with Gaussian noise
# Parameters: intercept b_0, gradient b_1, quadratic term b_2; Gaussian noise sigma, ysig.
# MCMC works on: theta=c(b_0, alpha=tan(b_1), b_2, log10(ysig)), a 1x4 vector.
# Prior PDFs:
# b_0:         N(mean=m, sd=s); m,s estimated from global properties of data
# alpha:       Uniform (0 to 2pi)
# b_2:         N(mean=m, sd=s); m,s estimated from global properties of data
# log10(ysig): Uniform (improper)

# Define covariance matrix of MCMC sampling PDF
sampleCov <- diag(c(0.1, 0.01, 0.01, 0.01)^2)
#lsfit <- lm(y ~ x + I(x^2), data=obsdata)
#summary(lsfit)
thetaInit <- c(-1, atan(-0.5), -1, log10(2.4))

# Run the MCMC to find postSamp, samples of the posterior PDF
set.seed(250)
allSamp <- metrop(func=logpost.quadraticmodel, thetaInit=thetaInit, Nburnin=2e4, Nsamp=2e5,
                   sampleCov=sampleCov, verbose=1e3, obsdata=obsdata)

##   1000 of  20000 +  2e+05  0.9169 
##   2000 of  20000 +  2e+05  0.8290 
##   3000 of  20000 +  2e+05  0.6824 
##   4000 of  20000 +  2e+05  0.5775 
##   5000 of  20000 +  2e+05  0.5201 
##   6000 of  20000 +  2e+05  0.4829 
##   7000 of  20000 +  2e+05  0.4597 
##   8000 of  20000 +  2e+05  0.4327 
##   9000 of  20000 +  2e+05  0.4112 
##  10000 of  20000 +  2e+05  0.4008 
##  11000 of  20000 +  2e+05  0.3951 
##  12000 of  20000 +  2e+05  0.3799 
##  13000 of  20000 +  2e+05  0.3751 
##  14000 of  20000 +  2e+05  0.3685 
##  15000 of  20000 +  2e+05  0.3580 
##  16000 of  20000 +  2e+05  0.3494 
##  17000 of  20000 +  2e+05  0.3417 
##  18000 of  20000 +  2e+05  0.3439 
##  19000 of  20000 +  2e+05  0.3418 
##  20000 of  20000 +  2e+05  0.3345 
##  21000 of  20000 +  2e+05  0.3333 
##  22000 of  20000 +  2e+05  0.3314 
##  23000 of  20000 +  2e+05  0.3301 
##  24000 of  20000 +  2e+05  0.3265 
##  25000 of  20000 +  2e+05  0.3233 
##  26000 of  20000 +  2e+05  0.3227 
##  27000 of  20000 +  2e+05  0.3219 
##  28000 of  20000 +  2e+05  0.3180 
##  29000 of  20000 +  2e+05  0.3175 
##  30000 of  20000 +  2e+05  0.3154 
##  31000 of  20000 +  2e+05  0.3191 
##  32000 of  20000 +  2e+05  0.3268 
##  33000 of  20000 +  2e+05  0.3307 
##  34000 of  20000 +  2e+05  0.3327 
##  35000 of  20000 +  2e+05  0.3367 
##  36000 of  20000 +  2e+05  0.3388 
##  37000 of  20000 +  2e+05  0.3378 
##  38000 of  20000 +  2e+05  0.3367 
##  39000 of  20000 +  2e+05  0.3348 
##  40000 of  20000 +  2e+05  0.3326 
##  41000 of  20000 +  2e+05  0.3323 
##  42000 of  20000 +  2e+05  0.3320 
##  43000 of  20000 +  2e+05  0.3358 
##  44000 of  20000 +  2e+05  0.3373 
##  45000 of  20000 +  2e+05  0.3363 
##  46000 of  20000 +  2e+05  0.3351 
##  47000 of  20000 +  2e+05  0.3351 
##  48000 of  20000 +  2e+05  0.3333 
##  49000 of  20000 +  2e+05  0.3322 
##  50000 of  20000 +  2e+05  0.3318 
##  51000 of  20000 +  2e+05  0.3324 
##  52000 of  20000 +  2e+05  0.3313 
##  53000 of  20000 +  2e+05  0.3292 
##  54000 of  20000 +  2e+05  0.3273 
##  55000 of  20000 +  2e+05  0.3256 
##  56000 of  20000 +  2e+05  0.3260 
##  57000 of  20000 +  2e+05  0.3269 
##  58000 of  20000 +  2e+05  0.3265 
##  59000 of  20000 +  2e+05  0.3247 
##  60000 of  20000 +  2e+05  0.3233 
##  61000 of  20000 +  2e+05  0.3233 
##  62000 of  20000 +  2e+05  0.3228 
##  63000 of  20000 +  2e+05  0.3217 
##  64000 of  20000 +  2e+05  0.3204 
##  65000 of  20000 +  2e+05  0.3205 
##  66000 of  20000 +  2e+05  0.3205 
##  67000 of  20000 +  2e+05  0.3198 
##  68000 of  20000 +  2e+05  0.3184 
##  69000 of  20000 +  2e+05  0.3199 
##  70000 of  20000 +  2e+05  0.3197 
##  71000 of  20000 +  2e+05  0.3189 
##  72000 of  20000 +  2e+05  0.3188 
##  73000 of  20000 +  2e+05  0.3178 
##  74000 of  20000 +  2e+05  0.3171 
##  75000 of  20000 +  2e+05  0.3160 
##  76000 of  20000 +  2e+05  0.3167 
##  77000 of  20000 +  2e+05  0.3157 
##  78000 of  20000 +  2e+05  0.3138 
##  79000 of  20000 +  2e+05  0.3123 
##  80000 of  20000 +  2e+05  0.3111 
##  81000 of  20000 +  2e+05  0.3104 
##  82000 of  20000 +  2e+05  0.3093 
##  83000 of  20000 +  2e+05  0.3090 
##  84000 of  20000 +  2e+05  0.3089 
##  85000 of  20000 +  2e+05  0.3088 
##  86000 of  20000 +  2e+05  0.3074 
##  87000 of  20000 +  2e+05  0.3063 
##  88000 of  20000 +  2e+05  0.3061 
##  89000 of  20000 +  2e+05  0.3054 
##  90000 of  20000 +  2e+05  0.3046 
##  91000 of  20000 +  2e+05  0.3062 
##  92000 of  20000 +  2e+05  0.3073 
##  93000 of  20000 +  2e+05  0.3083 
##  94000 of  20000 +  2e+05  0.3088 
##  95000 of  20000 +  2e+05  0.3091 
##  96000 of  20000 +  2e+05  0.3088 
##  97000 of  20000 +  2e+05  0.3092 
##  98000 of  20000 +  2e+05  0.3096 
##  99000 of  20000 +  2e+05  0.3090 
## 100000 of  20000 +  2e+05  0.3098 
## 101000 of  20000 +  2e+05  0.3094 
## 102000 of  20000 +  2e+05  0.3086 
## 103000 of  20000 +  2e+05  0.3085 
## 104000 of  20000 +  2e+05  0.3083 
## 105000 of  20000 +  2e+05  0.3085 
## 106000 of  20000 +  2e+05  0.3102 
## 107000 of  20000 +  2e+05  0.3103 
## 108000 of  20000 +  2e+05  0.3098 
## 109000 of  20000 +  2e+05  0.3094 
## 110000 of  20000 +  2e+05  0.3095 
## 111000 of  20000 +  2e+05  0.3094 
## 112000 of  20000 +  2e+05  0.3092 
## 113000 of  20000 +  2e+05  0.3087 
## 114000 of  20000 +  2e+05  0.3083 
## 115000 of  20000 +  2e+05  0.3076 
## 116000 of  20000 +  2e+05  0.3073 
## 117000 of  20000 +  2e+05  0.3070 
## 118000 of  20000 +  2e+05  0.3071 
## 119000 of  20000 +  2e+05  0.3063 
## 120000 of  20000 +  2e+05  0.3055 
## 121000 of  20000 +  2e+05  0.3042 
## 122000 of  20000 +  2e+05  0.3032 
## 123000 of  20000 +  2e+05  0.3023 
## 124000 of  20000 +  2e+05  0.3016 
## 125000 of  20000 +  2e+05  0.3013 
## 126000 of  20000 +  2e+05  0.3021 
## 127000 of  20000 +  2e+05  0.3023 
## 128000 of  20000 +  2e+05  0.3020 
## 129000 of  20000 +  2e+05  0.3020 
## 130000 of  20000 +  2e+05  0.3017 
## 131000 of  20000 +  2e+05  0.3014 
## 132000 of  20000 +  2e+05  0.3010 
## 133000 of  20000 +  2e+05  0.3011 
## 134000 of  20000 +  2e+05  0.3020 
## 135000 of  20000 +  2e+05  0.3030 
## 136000 of  20000 +  2e+05  0.3026 
## 137000 of  20000 +  2e+05  0.3028 
## 138000 of  20000 +  2e+05  0.3027 
## 139000 of  20000 +  2e+05  0.3032 
## 140000 of  20000 +  2e+05  0.3031 
## 141000 of  20000 +  2e+05  0.3027 
## 142000 of  20000 +  2e+05  0.3022 
## 143000 of  20000 +  2e+05  0.3020 
## 144000 of  20000 +  2e+05  0.3020 
## 145000 of  20000 +  2e+05  0.3017 
## 146000 of  20000 +  2e+05  0.3013 
## 147000 of  20000 +  2e+05  0.3010 
## 148000 of  20000 +  2e+05  0.3009 
## 149000 of  20000 +  2e+05  0.3017 
## 150000 of  20000 +  2e+05  0.3028 
## 151000 of  20000 +  2e+05  0.3041 
## 152000 of  20000 +  2e+05  0.3047 
## 153000 of  20000 +  2e+05  0.3051 
## 154000 of  20000 +  2e+05  0.3048 
## 155000 of  20000 +  2e+05  0.3040 
## 156000 of  20000 +  2e+05  0.3037 
## 157000 of  20000 +  2e+05  0.3030 
## 158000 of  20000 +  2e+05  0.3023 
## 159000 of  20000 +  2e+05  0.3021 
## 160000 of  20000 +  2e+05  0.3026 
## 161000 of  20000 +  2e+05  0.3030 
## 162000 of  20000 +  2e+05  0.3032 
## 163000 of  20000 +  2e+05  0.3026 
## 164000 of  20000 +  2e+05  0.3028 
## 165000 of  20000 +  2e+05  0.3031 
## 166000 of  20000 +  2e+05  0.3038 
## 167000 of  20000 +  2e+05  0.3041 
## 168000 of  20000 +  2e+05  0.3038 
## 169000 of  20000 +  2e+05  0.3030 
## 170000 of  20000 +  2e+05  0.3026 
## 171000 of  20000 +  2e+05  0.3021 
## 172000 of  20000 +  2e+05  0.3018 
## 173000 of  20000 +  2e+05  0.3015 
## 174000 of  20000 +  2e+05  0.3014 
## 175000 of  20000 +  2e+05  0.3015 
## 176000 of  20000 +  2e+05  0.3009 
## 177000 of  20000 +  2e+05  0.3006 
## 178000 of  20000 +  2e+05  0.3001 
## 179000 of  20000 +  2e+05  0.2997 
## 180000 of  20000 +  2e+05  0.2996 
## 181000 of  20000 +  2e+05  0.2998 
## 182000 of  20000 +  2e+05  0.3004 
## 183000 of  20000 +  2e+05  0.3002 
## 184000 of  20000 +  2e+05  0.2999 
## 185000 of  20000 +  2e+05  0.2996 
## 186000 of  20000 +  2e+05  0.2994 
## 187000 of  20000 +  2e+05  0.2992 
## 188000 of  20000 +  2e+05  0.2986 
## 189000 of  20000 +  2e+05  0.2981 
## 190000 of  20000 +  2e+05  0.2983 
## 191000 of  20000 +  2e+05  0.2995 
## 192000 of  20000 +  2e+05  0.2996 
## 193000 of  20000 +  2e+05  0.3000 
## 194000 of  20000 +  2e+05  0.2996 
## 195000 of  20000 +  2e+05  0.2998 
## 196000 of  20000 +  2e+05  0.3004 
## 197000 of  20000 +  2e+05  0.3011 
## 198000 of  20000 +  2e+05  0.3017 
## 199000 of  20000 +  2e+05  0.3024 
## 200000 of  20000 +  2e+05  0.3026 
## 201000 of  20000 +  2e+05  0.3039 
## 202000 of  20000 +  2e+05  0.3041 
## 203000 of  20000 +  2e+05  0.3044 
## 204000 of  20000 +  2e+05  0.3044 
## 205000 of  20000 +  2e+05  0.3043 
## 206000 of  20000 +  2e+05  0.3041 
## 207000 of  20000 +  2e+05  0.3038 
## 208000 of  20000 +  2e+05  0.3040 
## 209000 of  20000 +  2e+05  0.3042 
## 210000 of  20000 +  2e+05  0.3045 
## 211000 of  20000 +  2e+05  0.3045 
## 212000 of  20000 +  2e+05  0.3045 
## 213000 of  20000 +  2e+05  0.3048 
## 214000 of  20000 +  2e+05  0.3050 
## 215000 of  20000 +  2e+05  0.3049 
## 216000 of  20000 +  2e+05  0.3053 
## 217000 of  20000 +  2e+05  0.3060 
## 218000 of  20000 +  2e+05  0.3062 
## 219000 of  20000 +  2e+05  0.3057 
## 220000 of  20000 +  2e+05  0.3051

# 10^(allSamp[,1]+allSamp[,2]) is the unnormalized posterior at each sample
thinSel  <- seq(from=1, to=nrow(allSamp), by=100) 
postSamp <- allSamp[thinSel,]

# Plot MCMC chains and use density estimation to plot 1D posterior PDFs from these.
# Note that we don't need to do any explicit marginalization to get the 1D PDFs.
par(mfrow=c(4,2), mar=c(3.0,3.5,0.5,0.5), oma=0.5*c(1,1,1,1), mgp=c(1.8,0.6,0), cex=0.9)
parnames <- c(expression(b[0]), expression(paste(alpha, " / rad")), expression(b[2]), 
              expression(paste(log, " ", sigma)))
for(j in 3:6) { # columns of postSamp
  plot(1:nrow(postSamp), postSamp[,j], 
       type="l", 
       xlab="iteration", 
       ylab=parnames[j-2])
  postDen <- density(postSamp[,j], n=2^10)
  plot(postDen$x, postDen$y, 
       type="l",
       lwd=1.5,
       yaxs="i", 
       ylim=1.05*c(0,max(postDen$y)),
       xlab=parnames[j-2], 
       ylab="density")
  abline(v=thetaTrue[j-2], lwd=1.5, lty=3)
}

# Plot all parameter samples in 2D
par(mfcol=c(3,3), mar=c(3.5,3.5,0.5,0.5), oma=c(0.1,0.1,0.1,0.5), mgp=c(2.0,0.8,0))
for(i in 1:3) {
  for(j in 2:4) {
    if(j<=i) {
        plot.new()
      } else {
        plot(postSamp[,i+2], postSamp[,j+2], 
             xlab=parnames[i], 
             ylab=parnames[j], 
             pch=".")
    }
  }
}

# Find MAP and mean solutions.
# MAP = Maximum A Posteriori, i.e. peak of posterior.
# MAP is not the peak in each 1D PDF, but the peak of the 4D PDF.
# mean is easy, because samples have been drawn from the (unnormalized) posterior.
posMAP    <- which.max(postSamp[,1]+postSamp[,2]) 
thetaMAP  <- postSamp[posMAP, 3:6]
thetaMean <- apply(postSamp[,3:6], 2, mean) # Monte Carlo integration

# Overplot MAP solution with original data
par(mfrow=c(1,1), mar=c(3.5,3.5,2,1), oma=0.1*c(2,1,1,1), mgp=c(2.0,0.8,0), cex=1.0)
plotCI(obsdata$x, obsdata$y, 
       xlab="x",
       ylab="y", 
       uiw=10^thetaMAP[4],
       gap=0,
       col="royalblue",
       main="Quadratic fit of the Zeeman data")
xsamp <- seq(from=-3, to=2, length.out=500)
ysamp <- cbind(1,xsamp,xsamp^2) %*% as.matrix(modMat)
ysamp <- cbind(1,xsamp,xsamp^2) %*% as.matrix(c(thetaMAP[1], tan(thetaMAP[2]), thetaMAP[3]))
lines(xsamp, drop(ysamp), 
      lwd=1.5,
      col="red3") # MAP model
legend("topright", c("Data", "Fit"), col=c("royalblue", "red3"), lwd=3)

As can be seen from our data, the students forgot to take measurements in the region \(x \in (2.44, 3.0)\).

2. Run a Markov Chain Monte Carlo to predict peak height measurements at \(x_1 = 2.8\) mm and \(x_2 = 2.6\) mm.

x2.8 <- (2.8 - mean(x.in))/sd(x.in)
y2.8 <- thetaMAP[1] + x2.8*tan(thetaMAP[2]) + (x2.8^2)*thetaMAP[3] 

x2.6 <- (2.6 - mean(x.in))/sd(x.in)
y2.6 <- thetaMAP[1] + x2.6*tan(thetaMAP[2]) + (x2.6^2)*thetaMAP[3]
  
plot(x, y, 
     col="royalblue",
     xlab="x", 
     ylab="y", 
     main="Zeeman data + predicted data")
points(x2.8, y2.8, col="red3")
points(x2.6, y2.6, col="red3")
lines(xsamp, drop(ysamp), 
      lwd=1.5,
      lty=3,
      col="slategray3") 
legend("topright", c("Data", "Predictions"), col=c("royalblue", "red3"), lwd=3)

Exercise 2

The number of British coal mine disasters has been recorded from 1851 to 1962. By looking at the data it seems that the number of incidents decreased towards the end of the sampling period. We model the data as follows:

• before some year, we call \(\tau\), the data follow a Poisson distribution, where the logarithm of the mean value, \(\log \mu_t = b_0\), while for later years, we can model it as \(\log \mu_t = b_0 + b_1\);
• the dependence can be modeled as follows: \(y_t \sim \text{Pois}(\mu_t)\), where \(\log \mu_t = b_0 + b_1\text{Step}(t - \tau)\);
• implement the model in jags, trying to infer the parameters \(b_0\), \(b_1\) and \(\tau\);
• the step function is implemented, in BUGS, as step(\(x\)) and returns 1 if \(x \ge 0\) and 0 otherwise;
• assign a uniform prior to \(b_0\), \(b_1\) and a uniform prior in the interval \((1,N)\), where \(N = 112\) is the number of years our data span on;
• finally, here is our data:

data <− NULL
data$D <- c(4,5,4,1,0,4,3,4,0,6,3,3,4,0,2,6,3,3,5,4,5,3,
            1,4,4,1,5,5,3,4,2,5,2,2,3,4,2,1,3,2,1,1,1,1,
            1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,0,1,1,1,0,1,
            0,1,0,0,0,2,1,0,0,0,1,1,0,2,2,3,1,1,2,1,1,1,
            1,2,4,2,0,0,0,1,4,0,0,0,1,0,0,0,0,0,1,0,0,1,0,0)
data$N <− 112

• before running jags, assign an initial value to the parameters as follows: \(b_0 = 0\), \(b_1 = 0\) and \(\tau = 50\);
• explore the features of the chains and try to understand the effects of the burn-in, and thinning;
• plot the posterior distributions of the parameters and extract their mean values, and 95.

library(rjags)

## Loading required package: coda

## Linked to JAGS 4.3.0

## Loaded modules: basemod,bugs

library(coda)

data<- NULL

data$D <- c(4,5,4,1,0,4,3,4,0,6,3,3,4,0,2,6,3,3,5,4,5,3,
            1,4,4,1,5,5,3,4,2,5,2,2,3,4,2,1,3,2,1,1,1,1,
            1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,0,1,1,1,0,1,
            0,1,0,0,0,2,1,0,0,0,1,1,0,2,2,3,1,1,2,1,1,1,
            1,2,4,2,0,0,0,1,4,0,0,0,1,0,0,0,0,0,1,0,0,1,0,0)

data$N <- 112

model <- "model.bug"
parameters <- NULL
parameters$b0 <- 0
parameters$b1 <- 0
parameters$tau <- 50
jm <- jags.model(file = model, data = data, inits = parameters)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 115
##    Total graph size: 901
## 
## Initializing model

chain <- coda.samples(jm , c("b0", "b1", "tau"), n.iter=10000)
summary(chain)

## 
## Iterations = 1001:11000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD  Naive SE Time-series SE
## b0   1.141 0.09386 0.0009386       0.001847
## b1  -1.263 0.15461 0.0015461       0.002761
## tau 39.723 2.20381 0.0220381       0.058771
## 
## 2. Quantiles for each variable:
## 
##       2.5%    25%    50%    75%   97.5%
## b0   0.951  1.078  1.142  1.205  1.3180
## b1  -1.565 -1.366 -1.264 -1.158 -0.9611
## tau 36.106 37.984 39.928 40.841 45.1876

plot(chain, col="red3")

chain.df <- as.data.frame(as.mcmc(chain)) 
cat(sprintf("\n Correlation matrix: \n")) 

## 
##  Correlation matrix:

print(cor(chain.df))

##             b0          b1         tau
## b0   1.0000000 -0.57266116 -0.25404400
## b1  -0.5726612  1.00000000 -0.01303294
## tau -0.2540440 -0.01303294  1.00000000

t <- c(5, 10, 20, 40, 80)

for(i in 1:length(t)) {
    cat("Thinning =", t[i])
    chain <- coda.samples(jm, c("b0", "b1", "tau"), n.iter = 10000, thin = t[i])
    chain.df <- as.data.frame(as.mcmc(chain))
    
    b0.chain<-as.mcmc(chain.df["b0"])
    b1.chain<-as.mcmc(chain.df["b1"])
    tau.chain<-as.mcmc(chain.df["tau"])
    
    my.lags<-seq(0, 30)
    y1.b0 <- autocorr(b0.chain, lags=my.lags)
    y1.b1 <- autocorr(b1.chain, lags=my.lags)
    y1.tau <- autocorr(tau.chain, lags=my.lags)
    
    par(mfrow=c(3, 1))
    
    plot(my.lags, y1.b0,
         col="red3",
         xlab="lag", 
         ylab="ACF", 
         main=paste("b0 - thinning = ", t[i]))
    
    plot(my.lags, y1.b1, 
         col="red3",
         xlab="lag", 
         ylab="ACF", 
         main=paste("b1 - thinning = ", t[i]))

    plot(my.lags, y1.tau, 
         col="red3",
         xlab="lag", 
         ylab="ACF", 
         main=paste("tau - thinning = ", t[i]))

    plot(chain, col="red3")
    
    print(summary(chain))
}

## Thinning = 5

## 
## Iterations = 11005:21000
## Thinning interval = 5 
## Number of chains = 1 
## Sample size per chain = 2000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean     SD Naive SE Time-series SE
## b0   1.138 0.0932 0.002084       0.002084
## b1  -1.258 0.1503 0.003361       0.003643
## tau 39.710 2.0836 0.046591       0.056963
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%  97.5%
## b0   0.9524  1.078  1.140  1.203  1.319
## b1  -1.5548 -1.358 -1.256 -1.153 -0.965
## tau 36.1758 37.989 40.004 40.846 43.952
## 
## Thinning = 10

## 
## Iterations = 21010:31000
## Thinning interval = 10 
## Number of chains = 1 
## Sample size per chain = 1000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD Naive SE Time-series SE
## b0   1.138 0.09304 0.002942       0.002942
## b1  -1.263 0.15384 0.004865       0.004865
## tau 39.615 2.11662 0.066934       0.073731
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%   50%    75%   97.5%
## b0   0.9627  1.075  1.14  1.201  1.3180
## b1  -1.5485 -1.370 -1.26 -1.158 -0.9721
## tau 36.0753 37.923 39.82 40.807 44.2531
## 
## Thinning = 20

## 
## Iterations = 31020:41000
## Thinning interval = 20 
## Number of chains = 1 
## Sample size per chain = 500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD Naive SE Time-series SE
## b0   1.136 0.09026 0.004036       0.004036
## b1  -1.260 0.15747 0.007042       0.007561
## tau 39.574 2.24478 0.100389       0.100389
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%   97.5%
## b0   0.9567  1.078  1.141  1.200  1.3004
## b1  -1.5535 -1.376 -1.256 -1.149 -0.9493
## tau 36.0231 37.750 39.794 40.792 43.9902
## 
## Thinning = 40

## 
## Iterations = 41040:51000
## Thinning interval = 40 
## Number of chains = 1 
## Sample size per chain = 250 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD Naive SE Time-series SE
## b0   1.134 0.09047 0.005722       0.005722
## b1  -1.257 0.16375 0.010356       0.009239
## tau 39.705 1.99448 0.126142       0.154787
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%   97.5%
## b0   0.9847  1.071  1.133  1.191  1.3161
## b1  -1.5930 -1.368 -1.259 -1.140 -0.9524
## tau 36.1157 38.551 40.073 40.800 43.1573
## 
## Thinning = 80

## 
## Iterations = 51080:61000
## Thinning interval = 80 
## Number of chains = 1 
## Sample size per chain = 125 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD Naive SE Time-series SE
## b0   1.149 0.08342 0.007461       0.007461
## b1  -1.273 0.13754 0.012302       0.012302
## tau 39.752 2.08913 0.186858       0.186858
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%  97.5%
## b0   0.9738  1.096  1.164  1.207  1.291
## b1  -1.5301 -1.342 -1.267 -1.176 -1.045
## tau 36.2448 37.961 39.863 40.801 44.600

burn.in=c(10, 100, 200, 500, 1000, 2000, 4000)

for(i in 1:length(burn.in)) {
    jm <- jags.model(model, data, parameters, n.adapt = burn.in[i])
    chain <- coda.samples(jm, c("b0", "b1", "tau"), n.iter = 10000)
    
    cat("Burn-in = ", burn.in[i])
    
    par(mfrow=c(1, 1), fin=c(2, 2))
    
    plot(chain, col = "red3")
    
    chain.df <- as.data.frame(as.mcmc(chain))

    par(mfrow=c(3, 1), fin=c(4, 4))
    
    plot(chain.df$b0, chain.df$b1, 
         xlab = "b0", 
         ylab = "b1", 
         main = paste("b1 vs b0 - Burn-in =", burn.in[i]), 
         col = "red3")
    
    plot(chain.df$b0, chain.df$tau,
         xlab = "b0", 
         ylab = "tau",
         main = paste("tau vs b0 - Burn-in =", burn.in[i]),
         col = "red3")
    
    plot(chain.df$b1, chain.df$tau,
         xlab = "b1", 
         ylab = "tau", 
         main = paste("tau vs b1 - Burn-in =", burn.in[i]),
         col = "red3")
    
    print(summary(chain))
}

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 115
##    Total graph size: 901
## 
## Initializing model

## Warning in jags.model(model, data, parameters, n.adapt = burn.in[i]): Adaptation
## incomplete

## NOTE: Stopping adaptation
## 
## 
## Burn-in =  10

## 
## Iterations = 11:10010
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD  Naive SE Time-series SE
## b0   1.141 0.09433 0.0009433       0.001477
## b1  -1.261 0.15465 0.0015465       0.002290
## tau 39.735 2.24782 0.0224782       0.069223
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%   97.5%
## b0   0.9534  1.077  1.141  1.205  1.3221
## b1  -1.5638 -1.363 -1.261 -1.157 -0.9559
## tau 36.0620 37.960 39.953 40.843 46.0565
## 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 115
##    Total graph size: 901
## 
## Initializing model
## 
## Burn-in =  100

## 
## Iterations = 101:10100
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean     SD Naive SE Time-series SE
## b0   1.143 0.0932 0.000932       0.001761
## b1  -1.267 0.1541 0.001541       0.002736
## tau 39.725 2.1340 0.021340       0.055362
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%   97.5%
## b0   0.9579  1.080  1.145  1.207  1.3206
## b1  -1.5786 -1.369 -1.265 -1.163 -0.9736
## tau 36.1140 38.108 39.948 40.832 44.5496
## 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 115
##    Total graph size: 901
## 
## Initializing model
## 
## Burn-in =  200

## 
## Iterations = 201:10200
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD  Naive SE Time-series SE
## b0   1.137 0.09148 0.0009148       0.001686
## b1  -1.260 0.15147 0.0015147       0.002634
## tau 39.801 2.13088 0.0213088       0.061688
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%   97.5%
## b0   0.9508  1.078  1.138  1.201  1.3116
## b1  -1.5600 -1.361 -1.259 -1.157 -0.9623
## tau 36.1046 38.426 40.024 40.879 44.7881
## 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 115
##    Total graph size: 901
## 
## Initializing model
## 
## Burn-in =  500

## 
## Iterations = 501:10500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD  Naive SE Time-series SE
## b0   1.136 0.09406 0.0009406       0.001820
## b1  -1.255 0.15689 0.0015689       0.002872
## tau 39.721 2.14469 0.0214469       0.053148
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75% 97.5%
## b0   0.9485  1.073  1.138  1.201  1.32
## b1  -1.5660 -1.359 -1.253 -1.147 -0.95
## tau 36.1137 37.999 39.942 40.832 44.78
## 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 115
##    Total graph size: 901
## 
## Initializing model
## 
## Burn-in =  1000

## 
## Iterations = 1001:11000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD  Naive SE Time-series SE
## b0   1.140 0.09426 0.0009426       0.001897
## b1  -1.262 0.15667 0.0015667       0.002959
## tau 39.674 2.13750 0.0213750       0.054888
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%   97.5%
## b0   0.9498  1.078  1.142  1.204  1.3218
## b1  -1.5708 -1.369 -1.259 -1.156 -0.9568
## tau 36.1472 37.912 39.897 40.823 44.5590
## 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 115
##    Total graph size: 901
## 
## Initializing model
## 
## Burn-in =  2000

## 
## Iterations = 2001:12000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD  Naive SE Time-series SE
## b0   1.142 0.09368 0.0009368       0.001729
## b1  -1.264 0.15530 0.0015530       0.002861
## tau 39.737 2.12044 0.0212044       0.049448
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%  97.5%
## b0   0.9576  1.078  1.142  1.206  1.321
## b1  -1.5738 -1.366 -1.262 -1.159 -0.963
## tau 36.0929 38.235 39.977 40.857 44.376
## 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 115
##    Total graph size: 901
## 
## Initializing model
## 
## Burn-in =  4000

## 
## Iterations = 4001:14000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean      SD  Naive SE Time-series SE
## b0   1.139 0.09393 0.0009393       0.001767
## b1  -1.259 0.15661 0.0015661       0.002698
## tau 39.717 2.19728 0.0219728       0.051796
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%   97.5%
## b0   0.9499  1.077  1.141  1.202  1.3183
## b1  -1.5722 -1.363 -1.258 -1.153 -0.9507
## tau 36.1095 37.971 39.919 40.847 45.1527

Using different thinnings and burn-in we don’t see any important change in the results. Thus in this case it is more convenient to use a low thinning and a low burn-in (in order to get more samples).